clear;
% 猪肝样本的物理属性
rho = 1050; % 猪肝密度 (kg/m^3)
c = 3500; % 比热容 (J/kg·K)
kappa = 0.5; % 热导率 (W/m·K)

% 电极的加热区域体积 (椭球体积)
a = 0.01; % 长轴 (m)
b = 0.008; % 短轴 (m)
c_depth = 0.03; % 深度 (m)
V = (4/3) * pi * a * b * c_depth; % 椭球体积 (m^3)
P = 30; % 功率 (W)

% 离散化网格
r_max = 100; % 最大r
z_max = 100; % 最大z
dr = 0.1; % r的步长
dz = 1; % z的步长
dt = 0.01; % 时间步长
t_max = 8; % 最大时间

% 网格尺寸
r = -r_max:dr:r_max;
z = -z_max:dz:z_max;
t = 0:dt:t_max;

% 初始化温度场
T = zeros(length(r), length(z), length(t));

% 初始条件
T(:,:,1) = 25; % 初始温度 25℃

% 空间和时间步长相关常数
alpha = (kappa * dt) / (2 * rho * c);

% 进行时间步长迭代
for n = 1:length(t)-1
    % r方向的Crank-Nicolson系数
    for i = 2:length(r)-1
        if i > 2
            A(i-1,i-1) = 1 + 2*alpha / dr^2; % 主对角线
            A(i-1,i) = -alpha / dr^2; % 上对角线
            A(i-1,i-2) = -alpha / dr^2; % 下对角线
        end
    end
    
    % z方向的Crank-Nicolson系数
    for j = 2:length(z)-1
        if j > 2
            B(j-1,j-1) = 1 + 2*alpha / dz^2; % 主对角线
            B(j-1,j) = -alpha / dz^2; % 上对角线
            B(j-1,j-2) = -alpha / dz^2; % 下对角线
        end
    end
    
    % 使用Crank-Nicolson方法进行时间推进
    for i = 2:length(r)-1
        for j = 2:length(z)-1
            % 计算当前温度场和下一个时间步的温度场
            T(i,j,n+1) = T(i,j,n) + alpha * (T(i+1,j,n) - 2*T(i,j,n) + T(i-1,j,n)) ...
                                + alpha * (T(i,j+1,n) - 2*T(i,j,n) + T(i,j-1,n));
        end
    end
end
% 绘制最终的温度场
% 使用 imagesc 绘制颜色映射图
time_step = round((t_max-1) / dt);  % 最后一个时间步
T_at_time = T(:,:,time_step);

figure;
imagesc(r, z, T_at_time);  % 绘制二维温度场
axis xy;  % 设置坐标轴方向正确
xlabel('r (m)');
ylabel('z (m)');
title(['Temperature Distribution at t = ', num2str(time_step * dt), ' s']);
colorbar;  % 显示颜色条
